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Abstract. This paper is concerned with developing accurate and efficient numerical methods 
for one-dimensional fully nonlinear second order elliptic and parabolic partial differential equations 
(PDEs). In the paper we present a general framework for constructing high order interior penalty dis- 
continuous Galerkin (IP-DG) methods for approximating viscosity solutions of these fully nonlinear 
PDEs. In order to capture discontinuities of the second order derivative u xx of the solution u, three 
independent functions pi , p2 and P3 are introduced to represent numerical derivatives using various 
one-sided limits. The proposed DG framework, which is based on a nonstandard mixed formulation 
of the underlying PDE, embeds a nonlinear problem into a mostly linear system of equations where 
the nonlinearity has been modified to include multiple values of the second order derivative u xx . 
The proposed framework extends a companion finite difference framework developed by the authors 
in [9] and allows for the approximation of fully nonlinear PDEs using high order polynomials and 
non-uniform meshes. In addition to the nonstandard mixed formulation setting, another main idea 
is to replace the fully nonlinear differential operator by a numerical operator, which is consistent 
with the differential operator and satisfies certain monotonicity (called g-monotonicity) properties. 
To ensure such a g-monotonicity, the crux of the construction is to introduce the numerical moment, 
which plays a critical role in the proposed DG framework. The g-monotonicity gives the DG methods 
the ability to select the mathematically "correct" solution (i.e., the viscosity solution) among all pos- 
sible solutions. Moreover, the g-monotonicity allows for the possible development of more efficient 
nonlinear solvers as the special nonlinearity of the algebraic systems can be explored to decouple 
the equations. This paper also presents and analyzes numerical results for several numerical test 
problems which are used to guage the accuracy and efficiency of the proposed DG methods. 
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1. Introduction. Fully nonlinear partial differential equations (PDEs) refer to 
a class nonlinear PDEs which is nonlinear in the highest order derivatives of the 
unknown functions in the equations. Due to their strong nonlinearity, this class of 
PDEs are most difficult to analyze analytically and to approximate numerically. In the 
mean time, fully nonlinear PDEs arise in many applications such as antenna design, 
astrophysics, differential geometry, fluid mechanics, image processing, meteorology, 
mesh generation, optimal control, optimal mass transport, etc [8 , which calls for the 
development of efficient and reliable numerical methods for solving their underlying 
fully nonlinear PDE problems. 

This is the second paper in a series [9] which is devoted to developing finite 
difference (FD) and discontinuous Galerkin (DG) methods for approximating viscosity 
solutions of the following general one-dimensional fully nonlinear second order elliptic 
and parabolic equations: 

(1.1) F (u xx ,u x ,u,x) = 0, x e fl := (a, b), 
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and 

(1.2) u t + F(u xx ,u x ,u,t,x)=0, {x,t)eQ T :=nx(0,T), 

which are complemented by appropriate boundary and initial conditions. The goal 
of this paper is to design and implement a class of interior penalty discontinuous 
Galerkin (IP-DG) methods which are based on a nonstandard mixed formulation; the 
proposed IP-DG methods are named mlP-DG methods. For the ease of presenting 
the main ideas and avoiding the technicalities, in this paper we confine our attention 
to the one dimensional fully nonlinear second order PDE problem. The generalization 
and extension to the high dimensional case of the mlP-DG methods of this paper will 
be presented in a forthcoming work. In fact, it will be seen later that even in the one 
dimensional case, the construction and analysis of the proposed mlP-DG methods are 
already quite complicated. 

It is well known [5] that the primary challenges for approximating viscosity so- 
lutions of fully nonlinear PDEs are caused by the very notion of viscosity solutions 
themselves (see section [2] for the definition) . Unlike the notion of weak solutions 
for linear and quasilinear PDEs, the notion of viscosity solutions by design is non- 
variational, and, in general, viscosity solutions do not satisfy the underlying PDEs in a 
tangible sense. The non- variational nature of viscosity solutions immediately prevents 
any attempt to directly and straightforwardly construct Galerkin-type (including DG) 
methods for approximating fully nonlinear PDEs; in other words, nonlinearity in the 
highest order derivatives of the unknown function does not allow one to perform in- 
tegration by parts to transfer one order of derivatives to test functions as often done 
with linear and quasilinear PDEs. Another big challenge for approximating viscosity 
solutions of fully nonlinear PDEs is caused by the conditional uniqueness of viscosity 
solutions, namely, viscosity solutions may only be unique in a restricted function class. 
Requiring numerical solutions to stay or approximately stay in the same function class 
often imposes a difficult constraint for designing numerical methods. Finally, we like to 
mention that as expected, solving the resulting strong nonlinear (algebraic) systems, 
regardless which discretization method is used, is another difficult issue encountered 
with numerical fully nonlinear PDEs. 

The mlP-DG methods proposed in this paper aim to approximate viscosity solu- 



tions of (1.1) and (1.2) which belong to iT 1 (SI) in the spatial variable. We note that 
such a viscosity still does not satisfy the underlying PDEs in a tangible sense. We 
also mention that in order to approximate viscosity solutions that do not have H 1 
regularity in the spatial variable, we refer the reader to a companion paper |10j in 
which we propose another class of more complicated mixed discontinuous Galerkin 
that incorporates a local discontinuous Galerkin (LDG) approach instead of the IP- 
DG approach. Such an alternate LDG approach is also more appropriate when a more 
accurate approximation for u x is desired. 

Several novel ideas are utilized to design the mlP-DG methods in this paper which 
are briefly described below. Since integration by parts, which is the necessary tool 



for constructing any DG method, cannot be performed on equation (1.1), the first 
key idea is to introduce the auxiliary variable p := u xx and rewrite the original fully 
nonlinear PDE as a system of PDEs: 

(1.3) F(p, u x ,u, x) = 0, 

(1.4) p-u xx = 0. 

Unfortunately, since u xx may not exist for a viscosity solution u € H 1 (Q), the the 
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above mixed formulation may not make sense. To overcome this difficulty, the second 
key idea is to replace p := u xx by three possible values of u xx , namely, the left and 
right limits, as well as their average. Thus, we have 



(1.5) 
(1.6) 
(1.7) 



pi(x) - u xx (x~) = 0, 
p 2 (x) - u xx (x a ) = 0, 
p 3 (x) - u xx (x + ) = 0, 



where u xx (x a ) can be thought of as the arithmetic average of u xx (x~) and u xx (x + ). 
We remark that (1.5) and (1.7) can be regarded as two "one-sided" Poisson problems 
in u, and (1 1 .61) can be thought of as the "regular" Poisson problem. To incorporate 



the multiple values of u xx , equation (1.3) must be modified because F is only defined 
for a single value function p. To this end, we need the third key idea of this paper, 



which is to replace (1.3) by 



(1.8) 



F(pi,P2,P3,u x ,u,x) = 0, 



where F, which is called a numerical operator, should be some well-chosen approxi- 
mation to F. 

Natural questions that now arise are what are the criterions for F, and how can 
such a numerical operator F be constucted? These are two immediate questions which 
must be addressed. To do so, we need the fourth key idea of this paper, which is to 
borrow and adapt the notion of the numerical operators from our previous work [S] 
where a general finite difference framework has been developed for fully nonlinear 
second order PDEs. In summary, the criterions for F include consistency and g- 
monotonicity (generalized monotonicity) , for which precise definitions can be found 
in section|2j It should be pointed out that in order to construct the desired numerical 
operator F, a fundamental idea used in is to introduce the concept of the numerical 
moment, which can be regarded as a direct numerical realization for the moment term 
in the vanishing moment methodology introduced in |llj (also see (SJ section 4], [12 ). 



Finally, we need to design a DG discretization for the mixed system (1.5)-(1.8) to 



accomplish the goal. The fifth key idea of this paper is to use different numerical fluxes 



in the formulations of IP-DG methods for the "one-side" Poisson problems (1.5) and 



(1.7) as well as for the "regular" Poisson problem (1.6). We remark that, to the best 



of our knowledge, this is one of a few scenarios in numerical PDEs where the flexibility 
and superiority (over other numerical methodologies) of the DG methodology makes 
a vital difference. 

The remainder of this paper is organized as follows. In section [2] we collect some 
preliminaries including the definition of viscosity solutions, the definitions of the con- 
sistency and g-monotonicity of numerical operators, and the concept of the numerical 
moment. In section [3] we present the detailed formulation of mlP-DG methods for 
fully nonlinear elliptic equation ( 1.1 ) following the outline described above. In section 



[4] we consider both explicit and implicit in time fully discrete mlP-DG methods for 
fully nonlinear parabolic equation (1.2) based on the method of lines approach. The 



forward and backward Euler time-stepping schemes combined with the spatial mlP- 
DG methods will be specifically formulated. In section [5] we present many numerical 
experiments for the proposed mlP-DG methods and their fully discrete counterparts 



for the parabolic equation (1.2 1. These numerical experiments verify the accuracy of 



the proposed mlP-DG methods and also demonstrate the efficiency of these methods. 
Finally, we complete the paper with a brief summary and some concluding remarks 
in section [6j 
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2. Preliminaries. For a bounded open domain C R d , let B(Q), USC(Q) 
and LSC(fl) denote, respectively, the spaces of bounded, upper semi-continuous, and 
lower semicontinuous functions on n. For any v £ B(Q), we define 

v* (x) :— limsupu(y) and v*{x) := liminf v(y). 



y- 



y-tx 



Then, v* £ USC(Q) and u* £ LSC(Q), and they are called the upper and lower 
semicontinuous envelopes of v, respectively. 

Given a bounded function F : S dxd x R d x R x O -> R, where S dxd denotes the 
set of d x d symmetric real matrices, the general second order fully nonlinear PDE 
takes the form 

(2.1) F (D 2 u, Vu,u,x) = inO. 

Note that here we have used the convention of writing the boundary condition as a 
discontinuity of the PDE (cf. 2., p.274]). 

The following two definitions can be found in [7J [31 [2] . 

Definition 2.1. Equation (2.1) is said to be elliptic if for all (q, X,x) £ R d x 



R x f2 there holds 

(2.2) F(A, q, A, x) < F(B, q, A, x) VA,B £ S dxd , A> B, 

where A > B means that A — B is a nonnegative definite matrix. We note that when 
F is differentiable, the ellipticity also can be defined by requiring that the matrix 1? 
is negative semi-definite (cf. [7J p. 441]). 

Definition 2.2. A function u £ B(fl) is called a viscosity subsolution (resp. 



supersolution) of (2.1) if, for all tp £ C 2 (f2), if u* — ip (resp. — ip) has a local 



maximum (resp. minimum) at Xq £ Q, then we have 

F*(D 2 (p(x ),Vip(x ),u*(x ),x ) < 
(resp. F* (D 2 ip(xo) , V</?(xq), {%o), %o) > 0)- The function u is said to be a vis- 



supersolution of (2.1) 



cosity solution of (2.1) if it is simultaneously a viscosity subsolution and a viscosity 



We remark that if F and u are continuous, then the upper and lower * indices can 



be removed in Definition 2.2 The definition of ellipticity implies that the differential 



operator F must be non-increasing in its first argument in order to be elliptic. It turns 



out that ellipticity provides a sufficient condition for equation (2.1) to fulfill a maxi- 
mum principle (cf. [HE])- It is clear from the above definition that viscosity solutions 
in general do not satisfy the underlying PDEs in a tangible sense, and the concept of 
viscosity solutions is nonvariational. Such a solution is not defined through integra- 
tion by parts against arbitrary test functions; hence, it does not satisfy an integral 
identity. As pointed out in section[l] the nonvariational nature of viscosity solutions is 
the main obstacle that prevents direct construction of Galerkin-type methods, which 
are based on variational formulations. 

The following definitions are adapted from [9] in the case d = 1. 

Definition 2.3. 

(i) A function F : R 6 — > R is called a numerical operator. 
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(ii) A numerical operator F is said to be consistent ( with the differential operator 
F ) if F satisfies 

(2.3) liminf F(pi,p 2 ,p 3 , q x , X x , £1) > F*(p, q, X, £), 

Pfc — J-p. fc= 1 ,2 ,3 
Ql — — > — 

(2.4) limsup F(pi,P2,P3,qi,Xi,£i) < F*(p,q,X,£), 

where F* and F* denote respectively the lower and the upper semi-continuous 
envelopes of F. 

(iii) j4 numerical operator F is said to be g-monotone if F(pi,p 2 ,p 3 , q, A, £) is 
monotone increasing in p\ and p 3 and monotone decreasing in p 2 , that is, 

We note that the above consistency and g-monotonicity play a critical role in the 
finite difference framework established in [9j. They also play an equally critical role 
in the mlP-DG methods of this paper. We also note that in practice the consistency 
is easy to fulfill and to verify, but the g-monotonicity is not. In order to ensure the g- 
monotonicity, one key idea of [9 is to introduce the concept of the numerical moment 
to help. The following are two examples of so-called Lax-Friedrichs-like numerical 
operators [9]: 

(2.5) F 1 (p 1 ,p 2 ,p 3 ,q,X,£,) := F(p 2 ,q,X,£) + Qi(pi - 2p 2 + p 3 ) , 

(2.6) F 2 ( Pl: p 2: p 3:q ,X,0 := + ^ + P3 , q, X, g) + a 2 ( Pl - 2p 2 + p 3 ) , 



where a\ and a 2 are undetermined positive constants and the last term in (2.5 1 
and (2.6) is called the numerical moment. It is trivial to verify that F\ and F 2 are 
consistent with F . To ensure F\ to be g-monotone, we need 

dF 

(2.7) a > 



du" 

assuming adequate regularity for the operator F. We remark that it is natural to re- 
quire that F\ is decreasing in p 2 because by the definition of ellipticity, F is decreasing 
in u" . 

3. Formulation of mlP-DG methods for elliptic problems. We first con- 
sider the elliptic problem (1.1) with boundary conditions 

(3.1) u(a) = u a and u(b) = 

for two given constants u a and Ub- 

Let {xj}'j =0 C SI be a mesh for f2 such that xq — a and x.j — b. Define Ij = 
(xj-i,Xj) and hj — xj — Xj-± for all j = 1,2,..., J, ho = hj-\-% = and h = 
maxKjxjAj. Let 7~h denote the collection of the intervals {Ij}j =1 which form a 
partition of the domain J7. We also introduce the broken i? 1 -space 



ier h 
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and the broken L 2 -inner product 



(v,w) Th 



/ vwdx Vv, w e iJ 1 (7h). 

.7 = 1 Jl i 



For a fixed integer r > 1, we define the standard DG finite clement space V h C 
H^T^CL^Th) by 



yh 



where V r {I) denotes the set of all polynomials on / with degree not exceeding r. We 
also introduce the following standard jump and average notations: 

[vh(Xj)] := v h (xj) - v h (xf) for j = 1, 2, • • ■ , J - 1, 

K(x )] := -«fcOEo)> : = v h (xj); 

{v h (xj)} := ^(v h (xj) + v h (xffj for j = 1,2, ••• 
W^o)} := Vh(x ), {v h (xj)} := v h (xj). 
It is trivial to verify the following so-called "magic formulas" : 

(3.2) [v(xj)w(xj)\ = v(xJ)[w{ Xj )} + [v(xj)]w(xj), 

(3.3) [v(xj)w(xj)\ = Mx^jiwixj)} + [vix^iwixj)} 



J - 1, 



(3.4) 



[v(Xj)w(Xj)] = V(xj)[w(x j )} + [v{Xj)]w{Xj ). 



Let 70^ > for i = 1 , 2, 3 denote interior penalty parameters. It will be clear later 
that to avoid redundancy of three equations for pi,p2 and P3, we need to require that 
702 > max {701, 702}- Define the interior penalty terms 



(3.5) 



J 0i (V, W)=^2 — — [V (Xj)] [w (Xj)] 



3=0 



3,3+1 



for i = 1, 2, 3, where 

hj.j+i — max.{hj, hj+i} 



for j =0,1, 2,..., J. 



We now are ready to formulate our DG discretizations for equations (1.5 > — (1.8). 
First, for (fully) nonlinear equation (1.8) we simply approximate it by its broken 
L 2 -projection into V h , namely, 



(3.6) 
where 



ao(u h ,pi h ,p 2h ,p 3h ; <j>oh) = V<j3 0/i € V h , 

Oo(w,Pl,P2,P3;0o) = [F(P1,P2,P3,U',U, •),<&) 



75, 



Next, we discretize the three linear equations (1.5)— ( 1.7). Notice that for given 
"sources" {pi}f =1 , (1.5)-(1.7) are three (different) Poisson equations for u. Thus, we 
can use the standard IP-DG formulation for the Laplacian operator to discretize these 
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equations. However, there is a crucial distinction for doing so on the three equations, 



that is, we use, respectively, "magic formulas" (3.2 1, (3.3), and (3.4) when we add the 
local integration by parts formula to handle the jump terms at the interior nodes. To 
realize the above strategy, we define the bilinear forms bi : iJ 1 (77t) x -ff 1 (7/ l ) —> M by 

(3.7) bi(v,w) := (v',w') Th +v'(a)w(a) - ev(a)w'(a) - v'(b)w(b) 

+ ev(b)w'{b) + J 0i (v,w) Vv,w e H 1 (%) , i = 1,2,3, 

where e € {—1,0,1} is often called the "symmetrization" parameter [17]. Note, bi 
is symmetric if e = — 1, nonsymmctric if e = 1, and incomplete if e = 0. Using the 



bilinear forms bi, we define the following DG discretizations of (1.5)-(1.7) 
(3.8) ai(u h ,Pih;<t>ih) = M<Pih), V&h &V h , i = 1,2,3, 

where 

ai{u,px)(t>i) = (pi,<Ai)n + h{u,<t>i 
a 2 {u,p 2 ;<j> 2 ) = (j>2,<t>2)n + h(u,<t>2 
a3(w,p 3 ;03) = (P3, h)n + b t (u,4> 3 

and 

fitti) = Ma) -eti(a))u a + (-™- + 

V"0,1 7 V "J,.7+1 7 

for i = 1,2,3. 

In summary, our mlP-DG methods for the fully nonlinear Dirichlet problem (1.1 1, 
(1.2), and (3.1) are defined as seeking (u h ,p lh7 p 2 h,P3h) € [V h ] 4 such that (3.6) and 
~§ hold. 

We conclude this section with a few remarks. 

Remark 3.1. (a) Looking backwards, (3.8) provides a proper interpretation for 



.7-1 

E< 

i=i 


u'(xj )[(/>i( x j)] ~ 


e [u{x j )](j)' 1 {x j )), 


.7-1 

E( 

i=i 


^{ u '( x j)}[<h( x j)] 


- e [ u ( x i)]{<t>2{Xj)}), 


.7-1 

E( 


'u\xf)[Mxi)] - 


e[u{ Xj )]<&(xfj), 



(3. 



each ofpih, P2h> and p^h for a given function u^. Each pih defines a discrete second 
order derivative of u^. The functions Pih, P2h> an d P3h should be very close to each 
other if u xx exists; however, their discrepancies are expected to be large if u xx does 



not exist. p\h, P2h, and p^h defined by (3.8) can be regarded as high order extensions 



of their lower order counterparts defined in J3j. 

(b) It is easy to check that the three equations defined by (3.8) are linearly inde- 
pendent provided that 702 > max {701, 702}- 

(c) The reason for r ^ can be explained as follows. When r = 0, the piece- 
wise constant basis functions have piecewise zero derivatives on the given mesh. After 



eliminating the jump terms containing derivatives in (3.8), it is clear that the ability 
for pi and P3 to carry information from the left and the right, respectively, is lost. 
Furthermore, i/701 = 702 = 7o3> then a± = a 2 = 03, which in turn implies that on a 
uniform mesh pn t = p 2 ^ = p 3 h, they all are equal to the centered difference approxi- 
mation for the second order derivative ofu^. As a result, the numerical moment term 
vanishes and we are left with the standard three-point finite difference approximation 
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for (1.1 1 and (3.1 1. On the other hand, when r > 1, the numerical operator maintains 



the directional interpretations for pi and p 3; allowing the numerical operator to take 
advantage of the numerical moment. 

(e) Notice that (3.6)-( |3~8| ) is a nonlinear system of equations, with the nonlin- 
earity only appearing in ag. Thus, a nonlinear solver is necessary in implementing 
the above scheme. In section [5| an iterative method is used with initial guess given 
by projecting the secant line resulting from the boundary conditions into V h . Since a 
good initial guess is essential for most nonlinear solvers to converge, another possi- 
bility is to first linearize the nonlinear operator and solve the resulting linear system 
first. However, we show in our numerical tests that the simple initial guess works 
well in many cases. We suspect that the g-monotonicity of F enlarges the domain of 
"good" initial values over which the iterative method converges. 



Formulation of fully discrete mlP-DG methods for parabolic prob- 

The goal of this section is to extend our mlP-DG methods to solving the 



4. 

lems. 

initial-boundary value problem for (1.2) using the method of lines. Let the initial 
condition be given by 

(4.1) u(Q,x) = uo(x) V.t 6 SI 
and the boundary conditions be given by 

(4.2) u(t,a)=u a (t) ) u{t,b)=u b (t) Vt€(0,T).. 

For the ease of the presentation, we shall only consider the forward and backward 
Euler methods, although higher order time-stepping methods can also be formulated. 
For a fixed integer M > and let At = jj be the time-step size. Define t n := nAt 



for n = 0, 1, • • • , M. Then, the forward Euler method (in operator form) for ( 1.2 ) is 
defined by seeking u n+1 : SI — > R such that 



(4.3) 



,n+l 



= u n -AtF(u2 x ,u2,u n ,t n ,-) 



o 



and the standard backward Euler method (in operator form) for (1.2) is defined by 
seeking u n+1 : SI — >• R such that 



(4.4) u n+1 =u n - AtF(u n +\u n x + \u n +\t n +\ ■) 



for n = 0, 1, . . . , M — 1, where 



U° = Uq 



SI. 



The compatibility condition between the initial and boundary data immediately im- 
plies that 

(4.5) u"(o) = u a (t n ), u n (b) = u b (t n ), n = l,2,--.,M. 

We next apply the mlP-DG framework developed in the previous section to equa- 



tions (4.3) and (4.4) for their spatial discretizations. We present these two cases 



separately below because they require different treatments and involve different tech- 
nicalities. 
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4.1. Forward Euler method. It turns out that the forward Euler method is 
tricky to formulate because the variables u^ +1 and P™^ 1 {j = 1,2,3) are not deter- 
mined simultaneously. Instead, they are constructed sequentially. For a given itJJ, we 
first construct p™ h for j = 1,2,3 using (3.8). We then define u^ +1 to be a modified 
L 2 -projection of the right-hand side of (4.3 1. To take care of the boundary condition, 
we choose to enforce the boundary condition for u^ +1 weakly in the definition of the 
modified L 2 -projection. 

Specifically, for any v € L 2 (fl), we recall that the standard L 2 -projection VhV G 
V of v is defined by 

(4.6) (V h v, cf> h ) Th = (v, cf> h ) Th V<f> h e V h . 

For any v E C°(0), we define a modified L 2 -projection VhV <E V h of v by 

(4.7) (V h v,(/> h ) Th + i(^i)(o)^(o)+Pfc«{i)^(6)) 

= (v, 4>h) Th + (y(a)Ma) + v{b)4> h {b)) V(fo € V h , 

and the corresponding modified L 2 -projection operator 7\ ■ L 2 (0) fl C° (0) — > V' 1 . In 
the above definition, the boundary condition (4.2) is weakly enforced via a penalty 
technique which is due to Nitsche [To] . 

For a given function v e H 1 (7h) which satisfies the boundary condition (4.2), 
we define its three discrete "one-side" second order derivatives Q^w, Qh v > Qh v e 
(j = 1,2,3) using (3.8) as follows: 



(4.8) 



(^<M<z) ~ e<t>' h (a))u a (t) + (_^_0 h (6) + e #,(&))u 6 (t) 

.7-1 

- + ^(^'(^-)[^(^)] - eHxjM^xjj) E V h , 



3=1 



(Q^^) n = (^"^(a) - e^(a))«a(<) + + e& (&))«*(*) 



j-i 



(4.9) 



- b 2 (v, <f> h ) + ^ (K(^-)}[^(^)] - eK^OMfo)}) e v* 

3=1 
.7-1 

(4.io) -b 3 (v,ci )h ) + Y / (v , (4)[M^)}~4v(x 3 )W h (x+)) v<t> h ev h , 



3=1 



and the corresponding operators Q h , : if 1 (7/i) — > V' 1 . 

With the help of the operators Ql, Q^, Qt> an d T^hi we now define our fully 
discrete forward Euler method for the initial-boundary value problem (1.2), (4.1), 
d4~2l as follows: for n = 0, 1, . . . , M - 1, 



(4.11) 
(4.12) 



,71 + 1 



V h [u n h - MF(Q^ul, Q a h ul, QtK,{u n hx }, K}, f' 



u° h =V h u , 
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where {u^ x (x)} and {u^(x)} denote the average values of v% x and u\ at x. We remark 
that the average operation is needed because the approximation is only defined 
piecewisely on the mesh Th- Moreover, a projection is necessary at each time step 
since F may not belong to V h and the prescribed Dirichlet boundary condition must 
be enforced at each time step. 



as 



4.2. Backward Euler method. We rewrite the backward Euler scheme (4.4 1 

=u n - 1 



(4.13) u n + AtF( 

for n = 1, • • • , M, where u° = u a 



Clearly, (4.13) and (1.2 1 have the same form. Thus, the spatial discretization 



of the backward Euler scheme (4.4 1 is a straightforward adaptation of the mlP-DG 
framework for elliptic PDEs developed in section [3j To the end, we define a new 
numerical operator G by 

(4.14) G(pi,p 2 ,P3,u,t,x) :=u(x,t) + AtF(pi,p 2 ,P3,u x ,u,t,x) V(i,t)eflr. 



Then, analogous to the formulation for ( |1.1[) a nd (3.1 ), the fully discrete backward Eu- 
ler mlP-DG methods for ( fOj ), pa] ), andp^j) is defined by seeking (u%,p% h ,p% h ,pg h ) £ 
(V h ) 4 such that, for n = l,2,...,M, 



(4.15) ao^^h'PihAA^Oh) = «-\<Mr* 

(4-16) a i (<>Pth;<t>ih)=9i(t n ,<f>ih) 

(4.17) u° h = V h u . 

where 



aeV", * = 1,2,3, 



ao(t,u,pi,p 2 ,P3;^o) = (G(pi,P2,P3,U,t,-),4>o 



ai[u,pi; 



52(u'(xj)[Mxj)} 

3=1 


-fK^j^lh ))> 








]-eK^-)]{0' 2 (^)}) 


J'=l 










-e[^iM(4)), 



and 



= (j^ <M«) -eti(a))u a (t) + (-^- Mb) + e#(6)Wt) 

for i = 1, 2, 3. That is, a nonhomogeneous fully nonlinear elliptic problem is solved at 
each time step. 

5. Numerical experiments. In this section, we present a series of numerical 
tests to demonstrate the utility of the proposed mlP-DG methods for fully nonlinear 



PDEs of the types (1.1) and (1.2). In all of our tests we shall use uniform spatial 



meshes as well as uniform temporal meshes for the dynamic problems. To solve the 
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resulting nonlinear algebraic systems, we use the Matlab built-in nonlinear solver 
fsolve for the job. For the elliptic problems we choose the initial guess as the linear 
interpolant of the boundary data u a and u^. For dynamic problems, we let u° = VhUq, 
Pih = "LwOOi P°2h =_{u° hxx (x)}, and pl h = u° hxx (x+). Also, the initial guess for 
will be provided by u\ , and the initial guesses for p™ h , p% h , and p^ h will be provided 
by Pih 1 : P^h 1 j an< i Psih ' respectively. For convenience, we set e = for all tests. We 
remark that similar results can be obtained when e 0, and the actual benefit of the 
symmetrization parameter is unclear in the context of nonlinear algebraic systems. 



The role of a and the numerical moment will be further explored in section 5.3 

For our numerical tests, errors will be measured in the L°° norm and the L 2 norm, 
where the errors are measured at the current time step for the dynamic problems. 
For the dynamic test problems, we shall see that the lower order time discretization 
dominates the approximation error for reasonable time step size At. For the elliptic 
test problems and for the dynamic test problems where the error is dominated by the 
spatial discretizations, it appears that the spatial error is of order 0(h ), where 

^ [ r + 1, for r odd, 
1 r, for r even. 

Furthermore, we observe that when using odd order elements, the schemes exhibit 
optimal rate of convergence in both norms. 

5.1. Elliptic test problems. We first present the results for three test problems 



of type (1.1 ). Both Monge- Ampere and Bellman types of equations will be tested. 



Test 1. Consider the stationary Monge- Ampere problem 

-u 2 xx + 1 = 0, < x < 1, 
u(0) = 0, u(l) = i. 
It is easy to check that this problem has exactly two classical solutions: 

_i_ x \ 2 /\ ^2 

u (x) = -x , u (x) = ~^~ x > 

where u + is convex and u~ is concave. Note that u + is the unique viscosity solution 
which we want our numerical schemes to converge to. In section |5.3| we shall give 
some insights about the selectiveness of our schemes. 

We approximate the given problem using the linear element (r — 1) to see how 
the approximation converges with respect to h when the solution is not in the approx- 
imation space. The numerical results are shown in Figure |5.1| The results for the 



quadratic element (r = 2) are presented in Figure 5.12 We note that the approxima- 
tions using r = 2 are almost exact for each mesh size. This is expected since u + £ V h 
when r = 2. 

Test 2. Consider the problem 

-u 3 xx + \u x \+ S(x) = 0, -2<x<2, 
u(-2) = sin(4), u{2) = -sin(4), 



where 



S(x) = [2sign(x) cos(x 2 ) — Ax 2 sin(x|a;|)] 3 — 2|a;cos(a; 2 )|. 
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Test with N = 80, r= 1, and a= 2 



0.2 0.4 0.6 0.8 



r 


Norm 


h = 1/10 


h = 1/20 


h = 1/40 


h = 1/80 


Error 


Error Order 


Error Order 


Error Order 


1 


L 2 

L°° 


2.9e-03 
3.8e-03 


7.3e-04 2.00 
9.4e-04 2.00 


1.8e-04 1.99 
2.4e-04 1.99 


4.7e-05 1.97 
6.1e-05 1.96 



Fig. 5.1. Test 1: e = 0, a = 2, 701 = 703 = 1, and 702 = 1.1. 



This problem has the exact (viscosity) solution u(x) — sin(a;|a;|). Notice that the 
equation is nonlinear in both u xx and u x , and the exact solution is not twice differen- 
tiable at x = 0. The numerical results are shown in Figure [572) As expected, we can 
see from the plot that the error appears largest around the point x = 0, and both the 
accuracy and order of convergence improve as the order of the element increases. 

Test 3. Consider the stationary Hamilton- Jacobi-Bcllman problem 

inf 1 -6u xx + 9 2 x 2 u x + -u + S(x) 1 = 0, 1.2 < x < 4, 
o<eo)<i [ x J 



where 



S(x) 



u(1.2) = 1.44 In 1.2, u(4) = 16 In 4, 

4 ln(a;) 2 + 12 ln(x) + 9 - 8a; 4 ln(a;) 2 - 4a; 4 ln(x) 
4a; 3 [2 ln(x) + 1] ' 



It can be shown that the exact (viscosity) solution of this problem is given by u(x) = 
x 2 In a;, which occurs when 6*{x) — 2 X ^^h^ x J^T\ - ^ e s °l ve this problem using various 
order elements and record the numerical results in Figure |5.3[ which shows that our 
mlP-DG methods also can handle the Bellman-type fully nonlinear PDEs very well. 

5.2. Parabolic test problems. We now implement the proposed fully discrete 
forward and backward Euler mlP-DG methods for approximating fully nonlinear 



parabolic equations of the form (1.2). While the above formulation makes no at- 
tempt to formally quantify a CFL condition for the forward Euler method, our test 
problems generally require At — 0(h 2 ) to ensure the stability. In fact, the constant 
for the CFL condition appears to decrease as the order of the element increases. Below 
we implement both the implicit and explicit methods for each test problem. How- 
ever, we make no attempt to classify and compare the efficiency of the two methods. 
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r 


Norm 


h = 1 


h = 1/2 


h = l/4 


h = 1/8 


Error 


Error Order 


Error Order 


Error Order 


1 


L' 2 

L°° 


8.1e-01 
1.0e+00 


2.4e-01 1.73 
2.3e-01 2.14 


8.0e-02 1.60 
7.8e-02 1.58 


2.8e-02 1.52 
2.7e-02 1.54 


2 


L 2 

L°° 


l.le+00 
8.1e-01 


2.9c-01 1.88 
2.4e-01 1.76 


4.2e-02 2.78 
4.5e-02 2.40 


2.9e-02 0.56 
1.8e-02 1.30 


3 


L 2 

L°° 


6.4e-01 
4.9e-01 


2.7e-02 4.55 
3.1e-02 3.99 


1.4c-03 4.33 
1.6c-03 4.32 


6.5e-05 4.38 
9.1e-05 4.09 


4 


L 2 

L°° 


5.6e-02 
4.9e-02 


3.2e-03 4.14 
3.0e-03 4.02 


2.4e-04 3.72 
2.6e-04 3.56 


1.7e-05 3.83 
1.6e-05 4.02 


5 


L 2 

L°° 


2.3e-02 
2.1e-02 


8.5e-04 4.79 
9.3e-04 4.49 


1.5c-05 5.82 
1.8e-05 5.67 


2.4e-07 5.96 
2.6e-07 6.11 



Fig. 5.2. Test 2: e = 0, a = 4, 701 = 703 = 2, and 702 = 2.5. 



Instead, we focus on testing and demonstrating the usability of both fully discrete 
schemes and their promising potentials. For explicit tests, we record the parameter 
Kt which serves as the scaling constant for the CFL condition, so we have At — K t h 2 . 
For implicit tests, we record computed solution with various time step At. 

Test 4. Let = (0, 1), u a (t) — t 4 ,u b — | + t 4 , and uq{x) = \x 2 . We consider 



the problem O), Eli and (|4~2| with 



F(u xx ,u x ,u,t,x) = -u xx u+ -x 2 +t 4 - At 



3 



It is easy to verify that this problem has a unique classical solution u(x,t) = 0.5 x 2 + 
t 4 + 1. Notice that the PDE has a product nonlinearity in the second order derivative. 
The numerical results for the fully discrete forward Euler method are presented in 
Figure |5.4| and the results for the backward Euler method are shown in Figure |5.5| 
We observe that the errors for the backward Euler method are dominated by the 
relatively small size of the time step when compared to the forward Euler method. 
For smaller time step sizes, the errors are similar. However, the backward Euler 
method appears unstable for n t > 0.01. 

We now consider the error for the approximation resulting from using Euler time 
stepping methods. Note that the solution u is a quadratic in space. Letting r = 2, 
we limit the approximation error almost entirely to the time discretization scheme. 
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Test with N = 32, r= 1,and o=4 Test with N = 32, r= 4, and d = 4 




r 


Norm 


h = 2.8/4 


h = 2.8/8 


h = 2.8/16 


h = 2.8/32 


Error 


Error Order 


Error Order 


Error Order 


1 


L 2 

L°° 


3.5e-01 
3.9e-01 


9.8e-02 1.83 
1.2e-01 1.70 


2.6c-02 1.93 
3.4e-02 1.81 


6.6c-03 1.97 
9.0e-03 1.91 


2 


L 2 

L°° 


9.1e-03 
9.9e-03 


1.9e-03 2.28 
1.7e-03 2.53 


4.2e-04 2.18 
3.6e-04 2.23 


9.6e-05 2.11 
8.2e-05 2.15 


3 


L 2 

L°° 


3.5e-04 
5.1e-04 


2.7e-05 3.69 
4.2e-05 3.61 


1.9e-06 3.85 
3.3e-06 3.69 


4.2e-07 2.14 
3.7e-07 3.15 


4 


L 2 

L°° 


2.5e-05 
3.3e-05 


1.4e-06 4.14 
1.5c-06 4.46 


7.7e-08 4.19 
7.6e-08 4.30 


8.5e-09 3.18 
1.3c-08 2.51 



FlG. 5.3. Test 3: e = ? a = 4, 701 = 703 = 2, and 702 = 2.5. 



Explicit Test with N = 32, r= 1, and o= 2 at time = 1 Explicit rest with N= 32, r= 2, and d= 2 at time = 1 




0.2 0.4 0.6 n.8 1 n n.2 0.4 Q.B 0.8 1 



r 


Norm 


h = l/4 


h = 1/8 


h = 1/16 


h = 1/32 


Error 


Error Order 


Error Order 


Error Order 


1 


L 2 

L°° 


5.7e-03 
7.9e-03 


1.4e-03 1.98 
2.0e-03 1.99 


3.7e-04 1.99 
5.0e-04 1.99 


9.2e-05 1.99 
1.3e-04 1.99 


2 


L 2 

L°° 


3.3e-05 
4.5e-05 


8.2e-06 2.00 
l.le-05 2.00 


2.1e-06 2.00 
2.8e-06 2.00 


5.1e-07 2.00 
7.1c-07 2.00 


3 


L 2 

L°° 


3.3e-05 
4.5e-05 


8.2e-06 2.00 
l.lc-05 2.00 


2.1e-06 2.00 
2.8e-06 2.00 


5.1e-07 2.00 
7.1c-07 2.00 



Fig. 5.4. Test 4- Computed solutions at T = 1. Kt = 0.002, e = 0, a = 2, 701 = 703 = 2, and 
= 2.5. Note, the scheme is unstable for r = 2, 3 when Kt = 0.01. 
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2.5 
2 4 
2.3 
2.2 
2.1 
2 



0.2 0.4 0.6 0.8 1 



r 


Norm 


h = 1/4 


h = 


1/8 


h = 


1/16 






Error 


Error 


Order 


Error 


Order 


1 


L' 2 


4.4e-03 


9.6e-04 


2.20 


1.8e-04 


2.40 




L°° 


9.4e-03 


2.4e-03 


2.00 


5.9e-04 


2.00 


2 


L 2 


2.6e-04 


2.6e-04 


-0.00 


2.6e-04 


-0.00 




L°° 


3.6e-04 


3.6e-04 


-0.00 


3.6e-04 


-0.00 


3 


L 2 


2.6e-04 


2.6e-04 


-0.00 


2.6e-04 


-0.00 




L°° 


3.6e-04 


3.6e-04 


-0.00 


3.6e-04 


-0.00 



Fig. 5.5. Test 4- Computed solution at T = 1. At = 0.001, e = 0, a = 2, 701 = 703 = 2, and 
702 = 2.5. 



r 


Norm 


K t = 0.008 


Kt = 0.004 


K t = 0.002 


K t = 0.001 


Error 


Error Order 


Error Order 


Error Order 


2 


L' 2 

L°° 


8.2e-06 
l.le-05 


4.1e-06 1.00 
5.7e-06 1.00 


2.1e-06 1.00 
2.8e-06 1.00 


1.0e-06 1.00 
1.4c-06 1.00 



Fig. 5.6. Test 4-' Computed solutions at time T = 1. h = 1/16, e = 0, a = 2, 701 = 703 = 1, 
and 702 = 1.1. 




In fact, setting t = and solving the stationary form of the PDE, we have 



u h\\L 2 ({Q,l)) 



1.6 x 10" 



and ||u — u^II^oo^q^)) 1=3 2.4 x 10 



using the elliptic solver with h — 1/4, a = 2, 701 = 703 = 1, 702 = 1-1, and initial guess 
given by the secant line for the boundary data. Then, approximating the problem for 
varying Ai, we have the results recorded in Figure [BT6] for the forward Euler method 



and in Figure 5.7 for the backward Euler method. We observe that the convergence 
rate in time appears to have order 1 as expected. 

Test 5. Let Q = (0,2), u a (t) = 1, u b = e 2(t+1 ', and u (x) = e x . We consider the 
problem (l~2j), and (|4~2l with 



F(u xx ,u x ,u,t,x) 



:ln(u xx + l) + S(x,t), 



and 



S(x, t) = e [t+1)x (x - (t + 1) + lfe {t+1)x + l)) . 



(t+l)x 



It is easy to verify that this problem has a unique classical solution u(x,t) = e 
Notice that this problem is nonlinear in both u xx and u x . Furthermore, the exact 



16 XIAOBING FENG AND THOMAS LEWIS 



r 


Norm 


At = 1/10 


At = 1/20 


At = 1/40 


At = 1/80 


Error 


Error Order 


Error Order 


Error Order 


2 


L 2 

L°° 


2.4e-02 
3.3e-02 


1.3e-02 0.93 
1.7c-02 0.93 


6.4e-03 0.96 
8.8e-03 0.96 


3.2e-03 0.98 
4.5e-03 0.98 



Fig. 5.7. Test 4- Computed solutions at T = 1. h = 1/4, e = 0, a = 2, 701 = 703 = 1, and 
702 = 1-1- 



Explicit Test witri N = 32, r= 1,and c=4attlme = 0,5 Explicit Test with N= 32, r= 5, and o= 4 attlmc = 0.5 




r 


Norm 


h = 1/2 


h = 1/4 


h = 1/8 


h = 1/16 


Error 


Error Order 


Error Order 


Error Order 


1 


L 2 

L°° 


5.0e-01 
8.2e-01 


3.6e-02 1.73 
2.8e-01 1.57 


1.2e-02 1.32 
1.0e-01 1.47 


3.6e-03 1.67 
3.1e-02 1.69 


2 


L 2 

L°° 


4.5e-02 
6.0e-02 


1.2e-02 1.89 
1.4e-02 2.11 


3.3e-03 1.87 
3.6e-03 1.96 


8.7e-04 1.93 
9.0e-04 1.98 


3 


L 2 

L°° 


1.5e-03 
2.7e-03 


2.8e-04 2.39 
3.5e-04 2.97 


7.1e-05 1.98 
7.6e-05 2.21 


1.8e-05 1.98 
1.8e-05 2.05 


4 


L 2 

L°° 


1.2e-03 
1.3e-03 


2.9e-04 2.06 
3.0e-04 2.13 


7.2e-05 2.02 
7.3e-05 2.02 


1.8e-05 2.01 
1.8e-05 2.01 


5 


L 2 

L°° 


1.2e-03 
1.2e-03 


2.9e-04 2.00 
2.9e-04 2.00 


7.2e-05 2.00 
7.3e-05 2.00 


1.8e-05 2.00 
1.8e-05 2.00 



Fig. 5.8. Test 5: Computed solutions at time T = 3.10. Kt = 0.0025, e = 0, a = A, 701 = 
703 = 2, and 702 = 2.5. The method is not stable for Kt = 0.005. 



solution u cannot be factored into the form u(x, t) — G(t) Y(x) for some functions G 
and Y. The numerical results for the fully discrete forward Euler method are recorded 
in Figure |5.8| and the results for the backward Euler method are given in Figure |5.9| 
The error appears to be dominated by the low order time discretization given the 
relatively large value for At in the backward Euler test. However, using a smaller At 
for the forward Euler test, we were able to achieve a higher order of accuracy. We 
remark that even for At = 0.005/i 2 , the forward Euler scheme is not stable for h = | 
and r = 1. 

Test 6. Let il = (0, 2tt), u a (t) = 0, ui, = 0, and Uq(x) = sin(x). We consider the 



problem (1.2), fi~lL and (4.2) with 



F(u xx ,u x ,u,t, x) = - ^ min^ ^ ^A e u xx - c(x,i) cos(t) sin(a;) - sin(t) sin(a;)|, 
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ImplicitTest with N = 1 e, r= l,and o= 4 at time = 0.5 Implicit Test with N = 1 5, r= 3, and □ = 4 at time = 0.5 




r 


Norm 


h = 1/2 


h = 1/4 


h = 1/8 


Error 


Error Order 


Error Order 


1 


L' 2 

L°° 


4.2c-01 
8.3e-01 


1.4e-01 1.54 
2.4e-01 1.77 


4.6e-02 1.66 
7.9e-02 1.63 


2 


L 2 

L°° 


7.3e-02 
9.6e-02 


1.6e-02 2.21 
1.8e-02 2.41 


3.0e-03 2.40 
3.2e-03 2.49 


3 


L 2 

L°° 


2.8e-03 
5.6e-03 


7.8e-04 1.82 
8.5e-04 2.71 


9.1e-04 -0.22 
9.2e-04 -0.11 



Fig. 5.9. Test 5: Computed solutions at time T = 0.5. At = 0.0005, e = 0, a = 4, 701 = 703 = 
2, and 702 = 2.5. 



where A\ = 1, A2 = 5, and 

[l, if0<£<5 and 0<x<7ror^<t<7r and 7r < x < 2ir, 
I otherwise. 

It is easy to check that this problem has a unique classical solution u{x, t) = cos(t) sin(a:). 
Notice that this problem has a finite dimensional control parameter set, and the op- 
timal control is given by 

*•(*,*) 4 1 ' if f'*) = 1 ' 
v ' [2, if c(x,t) = 2. 

The numerical results are recorded in Figure |5.10| for the fully discrete forward Euler 
method and in Figure |5.11| for the backward Euler method. We observe that the 
accuracy of the implicit method appears to suffer from the lower order accuracy of 
the Euler method. For h = ? , the explicit method requires At 3.1 x 10 -4 , while 
the implicit method only needs At = 0.062. When At increases, the explicit method 
demonstrates instability. 

5.3. The role of the numerical moment. We now discuss the role and utility 
of the numerical moment in forming an appropriate numerical operator. Consider the 
stationary Monge- Ampere problem from Test 1, which has the following two solutions: 

_i_ x \ 2 /\ ^2 

u (x) = —x , u (x) — — x + x, 
w 2 ' v ' 2 

where u + is convex and u~ is concave. The solution u + is the unique viscosity solution. 
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Explicit Test with N = 32, r= 1,and 0= 2 at time = 3.1 Explicit Test with N = 32, r= 5, and □= 2 at time = 3.1 




r 


Norm 


h = ir/2 


h = tt/4 


h = ir/8 


h = tt/16 


Error 


Error Order 


Error Order 


Error Order 


1 


L 2 

L°° 


2.2e-01 
1.7e-01 


5.3e-02 2.07 
4.8e-02 1.87 


1.3c-02 2.02 
1.2c-02 1.98 


3.3e-03 2.01 
3.1e-03 1.99 


2 


L 2 

L°° 


6.0e-02 
6.4e-02 


1.6e-02 1.90 
1.5e-02 2.07 


4.2e-03 1.94 
3.5e-03 2.13 


l.le-03 1.97 
8.2e-04 2.09 


3 


L 2 

L°° 


7.4e-03 
8.0e-03 


6.9e-04 3.43 
5.6e-04 3.82 


1.4e-04 2.32 
1.0e-04 2.46 


3.5e-05 2.00 
2.3e-05 2.14 


4 


L 2 

L°° 


2.5e-03 
1.4e-03 


5.7e-04 2.10 
3.5e-04 2.01 


1.4e-04 2.03 
8.9e-05 1.98 


3.5e-05 2.01 
2.2c-05 1.99 


5 


L 2 

L°° 


2.2e-03 
1.4e-03 


5.6e-04 2.00 
3.6e-04 1.99 


1.4e-04 2.00 
8.9e-05 2.00 


3.5e-05 2.00 
2.2e-05 2.00 



Fig. 5.10. Test 6: Computed solutions at time T = 3.10. Kt = 0.002, e = 0, a = 2, 701 = 
703 = 2, and 702 = 2.5. 



To demonstrate the role of the numerical moment, we approximate the given 
problem using a > 0, a — 0, and a < 0. Notice that multiplying the PDE by 
— 1, we see that u~ is the unique viscosity solution of the equation with the operator 
F(u) = u 2 xx -l. Then, for a > 0, our scheme should converge to u + , and for a < our 
scheme should converge to u~ provided |a| is sufficiently large. However, for a = 0, 
the scheme may converge to either u + or vT depending on the initial guess used for 
the nonlinear solver. Note that while we cannot globally bound d U "F for the operator 
F(u") = 1 — (u") 2 , we can locally bound d U "F. Thus, the necessary magnitude for a 
to allow selective convergence depends on the initial guess and the solver. Without a 
global bound on d u »F, the numerical operator is only locally monotone. 

Let u be the linear interpolant of the boundary data and let the initial guess for 
u be given by — + | u~ and the initial guesses for pi be given by = 
for i — 1,2,3. Thus, the initial guess is closer to vT . From Figure 5.12 we see that 
the scheme converges to u + for a = 4 and the scheme converges to u~ for a = and 
a = —4 for the given parameters. If we change the initial guess to = |u + | u + , 
the scheme converges to u + for a — and a = 4 and the scheme converges to u~ for 
a = —4 for the given parameters. Furthermore, for u*- - 1 = u, fsolve does not find a 
root for a = 0, whereas the scheme converges for a — ±4. 

Therefore, the numerical moment plays two major roles. It allows the scheme to 
converge for a wider range of initial guesses, and it enables the scheme to address 
the issue of the conditional uniqueness of viscosity solutions. Given the form of the 
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Norm 


h = tt/2 


h = tt/4 


h = ir/8 


h = tt/16 


Error 


Error Order 


Error Order 


Error Order 


1 


L 2 

L°° 


1.7e-01 
1.5c-01 


4.9e-02 1.82 
4.4e-02 1.78 


1.4e-02 1.84 
1.3c-02 1.82 


4.8e-03 1.50 
4.1e-03 1.60 


2 


L 2 

L°° 


8.0e-02 
7.0e-02 


2.0e-02 2.00 
1.6e-02 2.14 


5.9e-03 1.76 
4.0e-03 1.98 


3.2e-03 0.87 
1.9c-03 1.06 


3 


L 2 

L°° 


l.le-02 
8.1e-03 


3.0e-03 1.91 
1.8e-03 2.16 


2.8e-03 0.09 
1.8e-03 0.01 


2.8e-03 0.00 
1.8e-03 0.00 



Fig. 5.11. Test 6: Computed solutions at time T = 3.10. At = 0.0062, e = 0, a = 2, 
701 = 703 = 2, and 702 = 2.5. 



numerical moment, a(pi — 2p2 + P3), these benefits are even more substantial given 
the way in which pi, P2, and ps are formed. The three variables only differ in their 
jump terms. When 701 = 702 = 703, the three different choices for the numerical 
fluxes (or jump terms) are all equivalent at the PDE level, and often the various jump 
formulations are presented as interchangeable when discretizing linear and quasilinear 
PDEs using the DG methodology. Yet, for our schemes for fully nonlinear PDEs, we 
see that the three different choices of the numerical fluxes all play an essential role 
at the numerical level when combined to form the numerical moment, even in the 
degenerative case where 701 = 702 = 703 which will be discussed below. 

The role of the numerical moment can heuristically be understood as follows when 
the numerical moment is rewritten in the form 

,2 (V\ ~ 2P2+P3 
<* h { p 

From here, we can see that the numerical moment acts as a centered difference approx- 
imation for u xxxx multiplied by a factor that tends to zero with rate O (h 2 ). Thus, 
at the PDE level, we are in essence approximating the nonlinear elliptic operator 

F {^XX ) ^X 7 ^5 ^) 

by the quasilinear fourth order operator F pi where 

Fp {u xxxx , 1L XX , 1L X , li, X) — p U xxxx -(- F {ti xx , U x , x) . 

In the limit as p — > 0, we heuristically expect the unique limit of the fourth or- 
der problem to converge to the unique viscosity solution of the second order prob- 
lem. Using a converging family of fourth order quasilinear PDEs to approximate a 
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Fig. 5.12. Test J: h = 1/10, e = 0, r = 2, 701 = 703 = 1.1, and 702 = 1.5. 



fully nonlinear second order PDE has previously been considered for PDEs such as 
the Monge-Ampere equation, the prescribed Gauss curvature equation, the infinity- 
Laplace equation, and linear second order equations of non-divergence form. The 
method is known as the vanishing moment method. We refer the reader to [121 [5] for 
a detailed exposition. 

In addition to the connection with the numerical moment to quasilinear fourth 
order PDEs, we also mention another benefit of the numerical moment. By the choice 
of a, we can enlarge the domain for which the numerical operator F is increasing in p\ 
and P3 and decreasing in p 2 ■ Since the definition of cllipticity is based on the mono- 
tonicity of the operator, and the issue of conditional uniqueness stems from whether 
the solution preserves the monotonicity of the operator, building monotonicity into 
the discretization is important when trying to preserve the nature of the operator we 
are approximating. 

We can demonstrate the power of the monotonicity of the numerical operator 
with two simple tests. For both tests, we shall again approximate the Monge-Ampere 
problem from Test 1 . However, now we let 701 = 702 = 703 • Then, we have p 2 = Pl ^ Pa , 
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Fig. 5.13. Plots of the norm of pi — 2p2 +P3 with a = 4, e = 0, and 701 = 702 = 703 = 2 at 
each iteration o/fsolve. 



which in turn implies that the equation for p 2 is redundant in the formulation and 
the numerical moment should be zero upon convergence to a root. 

For the first test, we again approximate the Monge- Ampere problem from Test 1 
while plotting the norm of p\ — 2p 2 +P3 after each iteration of f solve. From Figure [5.13| 
we can see that even though we expect the moment to be zero based on the redundancy 
of the equation for p 2 given the equations for p\ and p^, the Newton-based solver 
f solve treats p±, p 2 , and P3 as independent variables when searching for a root. The 
monotonicity of each variable appears to aid f solve in the search for a root. 

For the second test, instead of using fsolve, a Newton-based solver, for solving the 
nonlinear system of equations, we use the following splitting algorithm: 

Algorithm 5.1. 

(1) Pick an initial guess for u, p\, and p^. 



(2) Solve equation (3.6) for p 2 . 

(3) Solve equation (3.8) for i = 2 for u. 

(4) Solve equation (3.8) for i = 1 for pi. 

(5) Solve equation (3.8) for i = 3 for p^. 

(6) Repeat Steps 2-5 until the change in p 2 is sufficiently small. 

We observe that only Step (2) involves the use of a nonlinear solver. Each of 
Steps (3)-(5) only requires solving a linear system with a constant matrix that can 
be pre-factored. Thus, the above solver fully decouples the entire system of equations 
and minimizes the number of unknowns in the nonlinear system. Because this paper is 
concerned mainly with the discretization of fully nonlinear PDEs, we do not make an 
effort to compare solvers. The simple solver presented here is meant to demonstrate 
a potential benefit of using the numerical moment to create monotone numerical 
operators. 



We use Algorithm 5.1 with fsolve to execute Step (2) for Test 1. Let the initial 
guesses be given by u = u~ and p\ = p 2 = pz = —0.99. For p 2 = —0.99, F is increasing 
while F is decreasing for a > 0.99. Since F (—0.99) > and F is decreasing for p 2 > 
— 1 when a > 1, we expect the splitting algorithm will move away from the concave 
root p 2 = — 1 . The numerical results are presented in Figure |5T4"1 We note that even 
with the initial guess close to u~ , the solver, with the aid of the numerical moment, 
converges to u + . Similarly, the solver converges to u + when pi = p 2 = P3 > —0.99 
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are used as initial guesses. For initial guesses P\ = Pi = Pz < —1-0, the solver does 
not converge. Thus, we see that even for the above simple solver, the monotonicity of 
F provided by the numerical moment allows the scheme to either selectively converge 
to u + or diverge and find no solution. Hence, we again see the benefit of including 
the numerical moment when tackling the issue of conditional uniqueness for viscosity 
solutions. 

6. Conclusion. In this paper we present a general framework for constructing 
high order interior penalty discontinuous Galerkin methods for approximating viscos- 
ity solutions of fully nonlinear second order elliptic and parabolic PDEs. The proposed 
framework extends the (second order) finite difference framework developed by the 
authors in [9] to a more flexible DG framework, allowing approximating fully nonlin- 
ear PDEs using high order polynomials and non-uniform meshes. Various numerical 
experiments are provided to show the performance of the proposed methodology. The 
proposed DG framework is based on a nonstandard mixed formulation of the under- 
lying fully nonlinear PDE. In order to capture discontinuities of the second order 
derivative u xx of the solution u, three independent functions p\,P2 and p% are intro- 
duced to accomplish the goal, where p\ and p^ measure the left and right limits of 
u". If u" is discontinuous, p\ and p 3 can be used to gain insight into the discontinuity 
upon convergence. Thus, the methodology has the ability to capture some of the 
more interesting aspects of the viscosity solutions. The proposed mlP-DG methodol- 
ogy takes the most important aspects of the companion finite difference framework of 
[9] and extends them in multiple directions. For example, by adopting and expand- 
ing the idea of numerical operators, the mlP-DG formulation allows for even more 
flexibility than finite difference methods in construction. 

The proposed mlP-DG discretizations touch the inner core and make use of the 
full potential of the DG methodology. This is because there is a natural match among 
the three choices of possible numerical fluxes and the three numerical second order 
derivatives pi, P2, and P3, and the flexibility of DG methods allows the implantation 
of this connection into the formulation of the proposed mlP-DG methods. 

Like in the finite difference framework, the g-monotonicity (generalized mono- 
tonicity) and the numerical moment play a central role in the proposed mlP-DG 
framework. The g-monotonicity gives the mlP-DG methods the ability to select the 
mathematically "correct" solution (i.e., the viscosity solution) among all possible solu- 
tions, and the numerical moment is the catalyst which facilitates the g-monotonicity of 
the proposed mlP-DG methods. Moreover, the g-monotonicity allows for the possible 
development of more efficient (than generic Newton) solvers. The special nonlinearity 
of the algebraic systems can be explored to decouple the equations as seen in Algo- 



rithm 5.1 We believe that one of the main strengths of the mlP-DG formulation 
presented in this paper lies in the way in which the discretization handles the non- 
linearity. The discretization takes a nonlinear problem and embeds it into a mostly 
linear system of equations where the nonlinearity has been modified to ensure g- 
monotonicity. The added monotonicity can theoretically enlarge the domain of valid 
initial guesses over which a solver will converge. Thus, the weak coupling with linear 
equations is only a small penalty for the added structure in the nonlinearity. 

We also remark that the role and benefit of the symmetrizing parameter is unclear 
for nonlinear systems of equations. When 70 is sufficiently large, we observe that 
numerical results seem independent of the choice for e. However, when the penalty 
constant 70 is not sufficiently large, the inclusion of e can be detrimental to the 
approximation. For small 70, the approximation is allowed to have larger jumps 
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6.3e-04 1.99 
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Fig. 5.14. Test 1 is solved using Algorithm \5 . l\ with r = 1, e = 0, and 701 = 702 = 703 = 2. For 
a > 1.1, the scheme converges to u + . For a < 0.99, the scheme converges to u~ . When a = 1.0, 
the scheme converges to u + for h > and the scheme converges to u~ for h < wfj. 



occur. When the jumps become too large, the effect from having e present becomes 
exaggerated, and the overall accuracy of the approximation begins to suffer beyond 
just the presence of jumps. For elliptic problems, we can see the formation of a 
boundary layer. For dynamic problems, we see the approximation actually diverging 
(almost instantaneously) along the interior of the domain. As expected, when 70 
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increases and becomes sufficiently large, these phenomena disappear. Thus, at a 
numerical level, the presence of the symmetrizing constant e seems important, even 
though at the continuous level for the PDE, the symmetrization terms all become 
zero. 

Conceptually, the mlP-DG framework presented in this paper can be easily ex- 
tended to the high dimensional fully nonlinear PDE problems, the detailed exposition 
will be given in a forthcoming paper. On the other hand, the proposed mlP-DG 
framework may not work in the case when the viscosity solution does not belong to 
In such a case, a more involved mixed local discontinuous Galerkin (mLDG) 
framework must be invoked. We refer the interested reader to [TU] for a detailed 
exposition. 
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